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ABSTRACT 

Boussinesq-type  equations  for  weakly  nonlinear,  weakly  dispersive  waves 
have  been  used  extensively  to  model  wave  shoaling  on  beaches.  Deterministic 
Boussinesq  models  cast  in  the  form  of  coupled  evolution  equations  for  the 
amplitudes  and  phases  of  discrete  Fourier  modes  (Freilich  and  Guza,  1984)  describe 
the  shoaling  process  accurately  for  arbitrary  incident  wave  conditions,  but  are 
numerically  cumbersome  for  predicting  the  shoaling  evolution  of  continuous  spectra 
of  natural  wind-generated  waves.  Here  an  alternative  stochastic  formulation  of  a 
Boussinesq  model  (Herbers  and  Burton,  1996,  based  on  the  closure  hypothesis  that 
phase  coupling  between  quartets  of  wave  components  is  weak)  is  implemented  that 
predicts  the  evolution  of  a  continuous  frequency  spectrum  and  bispectrum  of  waves 
normally  incident  on  a  gently  sloping  beach  with  straight  and  parallel  depth 
contours.  The  general  characteristics  of  the  model  are  examined  with  numerical 
simulations  for  a  wide  range  of  incident  wave  conditions  and  bottom  profiles. 
Stochastic  and  deterministic  Boussinesq  model  predictions  are  compared  to  field 
observations  from  a  cross-shore  transect  of  bottom  pressure  sensors  deployed  on 
a  barred  beach  near  Duck,  NC,  during  the  recent  DUCK94  Experiment.  Predictions 
of  the  two  models  are  similar  and  describe  accurately  the  observed  nonlinear 
shoaling  transformation  of  wave  spectra. 
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I.   INTRODUCTION 

Wind-generated  surface  gravity  waves  are  the  principal  driving  force  of  nearshore 
fluid  motions  (e.g.,  longshore  currents,  rip  currents,  and  undertow)  and  sediment  transport 
(e.g.,  erosion  and  accretion  of  beaches,  and  the  formation  of  bars  and  cusps).  As  waves 
shoal  onto  beaches,  amplitudes  increase,  wavelengths  decrease,  and  propagation  directions 
refract  towards  normal  incidence.  These  linear  propagation  effects  are  readily  observed 
and  well  understood.  Additionally,  pronounced  nonlinear  effects  in  shallow  water  cause 
a  dramatic  transformation  of  wave  shapes  from  initially  smooth,  nearly  sinusoidal  profiles, 
to  asymmetric,  pitched  forward  profiles  characteristic  of  near-breaking  waves.  The 
mechanism  for  this  transformation  is  nonlinear  triad  interactions  in  which  two  primary 
wave  components  with  frequencies  o)1  and  (D2  excite  a  secondary  wave  component  with 
the  sum  (o> 2+o>2)  or  difference  (<yrto2)  frequency.  The  nonlinear  energy  transfers  in 
these  interactions  not  only  broaden  the  spectrum  but  also  change  the  statistical  properties 
of  the  waves.  Whereas  the  primary  wave  components  incident  from  deep  water  are 
approximately  statistically  independent,  the  newly  formed  secondary  components  are 
"phase  locked"  to  the  primary  waves  that  excite  them.  Even  relatively  weak  secondary 
components  significantly  change  the  shapes  of  waves  in  shallow  water.  While  the  incident 
waves  and  the  nonlinearly  excited  higher  frequency  waves  are  predominantly  dissipated 
in  the  surf  zone,  the  nonlinearly  excited  lower  frequency  (infragravity)  wave  components 
reflect  from  the  beach  and  often  dominate  wave  runup  at  the  shoreline.  Accurate 
prediction  of  the  nonlinear  shoaling  transformation  of  waves  is  critical  to  both  naval  (e.g., 
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amphibious  landings  and  mine  warfare)  and  civilian  (e.g.,  mitigation  of  beach  erosion) 
operations  in  the  littoral  zone. 

In  deep  water  (kJi  »  1,  where  k  is  the  wavenumber  and  h  is  the  water  depth)  and 
intermediate  depths  (kJi  =  0(1))  triad  interactions  are  non-resonant.  The  nonlinearly 
excited  secondary  waves  remain  small  ("bound"  waves)  and  are  well  described  by  finite 
depth  theory  based  on  Stokes  perturbation  expansion  for  small  wave  steepness  (Phillips 
1960;  Hasselmann  1962;  Herbers  et  al.,  1992,1994;  and  many  others).  In  shallow  water 
(kJi  <  1),  where  triad  interactions  are  near-resonant,  finite  depth  theory  is  valid  only  for 
small  values  of  the  Ursell  number,  Ur  =  a/Kh  (Ursell,  1953)  where  a  is  the  wave 
amplitude,  and  this  condition  is  usually  violated  on  natural  beaches. 

Models  for  waves  in  shallow  water  are  usually  based  on  the  Boussinesq  equations 
which  assume  that  a/h  (nonlinearity)  and  {xh)  (dispersion)  are  both  small  and  of  the 
same  order  (i.e.,  Ur  =0(1)).  These  equations  are  surprisingly  robust,  even  for  large  Ur 
values  typically  observed  in  near-breaking  waves.  Peregrine  (1967)  extended  the 
Boussinesq  equations  to  varying  depth,  and  these  equations  form  the  basis  of  most  wave 
shoaling  models.  Freilich  and  Guza  (1984)  developed  a  frequency  domain  (i.e.,  neglecting 
directional  spreading  effects)  Boussinesq  model  that  predicts  accurately  the  energy 
transfers  to  higher  frequencies  and  associated  wave  shape  changes  on  natural  beaches 
(Elgar  and  Guza,  1985a;  Elgar  et  al.,  1990a).  This  deterministic  model  is  initialized  at  an 
offshore  boundary  by  a  discrete  Fourier  representation  of  the  incident  waves.  A  coupled 
set  of  evolution  equations  for  the  amplitudes  and  phases  of  the  Fourier  modes  is  solved 
numerically  to  evaluate  the  shoaling  transformation  of  the  wave  train.  The  results  of  these 


integrations  are  subsequently  averaged  over  many  realizations  with  random  initial 
amplitudes  and  phases  to  obtain  spectral  statistics.  Although  the  shoaling  evolution  of 
wave  spectra  is  accurately  reproduced,  this  method  is  cumbersome  for  practical 
applications,  requiring  large  computing  resources  and  a  detailed  specification  of  incident 
wave  conditions  at  the  offshore  boundary  that  is  often  not  available.  In  addition,  the 
extension  of  this  approach  to  include  directional  spreading  effects  (i.e.,  two  dimensions) 
is  far  from  straightforward  owing  to  the  large  number  of  modes  required  in  discrete  mode 
simulations  of  continuous  frequency-directional  wave  spectra. 

The  initialization  of  wave  shoaling  models  requires  local  offshore  measurements 
(e.g.,  wave  following  buoys)  or  predictions  from  a  global  or  regional  wind-wave 
generation  model.  These  larger  scale  models  (e.g.,  WAM;  The  WAMDI  group,  1988)  are 
typically  stochastic  formulations  based  on  an  energy  balance,  and  predict  wave  spectra 
rather  than  individual  wave  profiles.  Routine  measurement  systems  provide  similar 
statistical  information. 

Recently,  Herbers  and  Burton  (1996)  derived  an  alternative  stochastic  formulation 
of  Boussinesq  shoaling  evolution  equations,  for  directionally  spread  waves  propagating 
over  a  gently  sloping  beach  with  straight  and  parallel  contours.  Based  on  the  closure 
hypothesis  that  phase-coupling  between  quartets  of  wave  components  is  weak,  Herbers 
and  Burton  derived  a  coupled  set  of  evolution  equations  for  the  wave  spectrum  and 
bispectrum.  The  bispectrum  describes  the  degree  of  coupling  and  the  phase-relationship 
in  triads  of  nonlinearly  interacting  wave  components  (Hasselmann  et.  al.,  1963),  and  is 
used  extensively  in  other  nonlinear  process  studies  (e.g.,  economics,  brain-wave  emission, 


and  plasma  physics).  In  deep  water  and  intermediate  depths,  the  bispectrum  is  completely 
determined  by  the  local  spectrum,  and  enables  the  detection  of  relatively  weak  phase- 
coupled,  forced  secondary  waves  that  are  concealed  in  the  spectrum  by  more  energetic 
freely  propagating  primary  waves  (e.g.,  Herbers  et  al.,  1992,  1994).  In  shallow  water,  the 
bispectrum  evolves  strongly  and  describes,  in  a  statistical  sense,  the  shape  evolution  of 
shoaling  waves  (e.g.,  Elgar  and  Guza,  1985b;  Elgar  et  al.,  1990a). 

Here,  a  numerical  implementation  of  a  one-dimensional  stochastic  Boussinesq 
model  is  presented  in  which  directional  spreading  effects  are  neglected.  This  formulation 
allows  for  simple  illustration  of  stochastic  model  characteristics  and  comparisons  to  field 
data  and  existing  one-dimensional  models.  Energy  transfers  to  higher  frequencies  in  sum 
triad  interactions  are  insensitive  to  directional  spreading  angles  of  incident  waves  (e.g., 
Herbers  and  Burton,  1996),  and  thus  can  be  predicted  accurately  with  a  one-dimensional 
model.  However,  energy  transfers  to  infragravity  frequencies  in  difference  triad 
interactions  are  significantly  reduced  for  large  directional  spreading  angles  (Herbers  and 
Burton,  1996).  Additionally,  the  reflection  from  shore  and  refractive  trapping  of 
infragravity  waves  (Herbers  et  al.,  1995)  is  neglected.  Hence,  infragravity  waves  are  only 
crudely  represented  in  the  present  model  formulation. 

Following  a  review  of  the  stochastic  formulation  of  Boussinesq  wave  shoaling 
equations,  the  numerical  model  implementation  is  described  in  Chapter  II.  The 
dependence  of  wave  shoaling  evolution  on  the  nonlinearity,  spectral  shape  of  incident 
waves  and  the  beach  profile  is  examined  through  numerical  simulations  in  Chapter  III. 
Stochastic  and  detenninistic  (the  Freilich  and  Guza  model)  Boussinesq  predictions  are 


compared  to  field  data  collected  on  a  natural  beach  near  Duck,  NC  in  Chapter  IV, 
followed  by  conclusions  in  Chapter  V. 


II.  A  STOCHASTIC  BOUSSINESQ  MODEL 

Herbers  and  Burton  (1996)  derived  a  stochastic  formulation  of  Boussinesq  wave 
shoaling  equations  for  directionally  spread  waves  propagating  over  a  beach  with  straight 
and  parallel  depth  contours.  Under  the  third-order  closure  hypothesis  that  phase-coupling 
between  quartets  of  wave  components  is  weak,  the  statistical  properties  of  the  waves  are 
described  by  a  coupled  set  of  evolution  equations  for  the  frequency  (co) -alongshore 
wavenumber  (/)  spectrum  E{o),t)  and  bispectrum  B((x)',l',G}-(*i',l-lr).  If  directional 
spreading  is  neglected  (i.e.,  /  =  0),  these  equations  (22a  and  22b  in  Herbers  and  Burton, 
1996)  reduce  to: 


d  1  dh  3g)       r     ,       i  /,      / 

—E((D)  = £(g>)  +  flM{B(w,u-to)}dG> 

dx  2hdx  2A3V/2J- 


(1) 


d         /  / 

— 5(0)  ,o)-(o  )  = 
dx 


3  dh       Ji     co  (g)-g>  )co 
~4hdx       l  ^3/2 


B((x>  ,0)-G)  ) 


(2) 


-  i- 


2A3V/2L 


co  £((o-co  )£(o))+(co-a)  )£(<o  )£(a))-a)£(o)  )£(g>-g)  ) 


where  the  x-axis  points  onshore,  E(a>)  and  B(d) ',  o)-co  0  are  the  frequency  spectrum  and 
bispectrum,  h(x)  is  the  water  depth,  g  is  gravity,  and  IM{  }  indicates  the  imaginary  part. 
The  integrals  of  the  (density)  spectrum  E  and  bispectrum  B  over  all  frequencies  yield  the 
mean  square  <ti  >  and  mean  cube  <r|  >  of  the  surface  elevation  r\(x,t): 


<ti2>  =    /*£(<*)  )d(*> 


<T|  >  =   f  [B((d  ,G)-G)  )doi  d(x) 


The  first  term  on  the  right-hand  side  of  equations  (1)  and  (2)  represents  linear 
shoaling  effects.  The  nonlinear  transfers  in  the  energy  spectrum  are  controlled  by  the 
imaginary  part  of  the  bispectrum  (the  integral  on  the  right-hand  side  of  equation  (1)).  The 
energy  product  terms  in  equation  (2)  represent  the  changes  in  the  imaginary  part  of  the 
bispectrum  due  to  the  three  possible  nonlinear  interactions  (one  sum  interaction  and  two 
difference  interactions)  within  the  (a)',co-co',(o)  triad.  The  second  term  on  the  right-hand 
side  of  equation  (2)  represents  the  detuning  of  the  interactions  from  resonance  (i.e., 
changes  in  the  phase  of  the  bispectrum  owing  to  weak  dispersion).  In  the  limit  of  small 
amplitudes  and  bottom  slope,  steady  solutions  for  E(co)  and  B(o)',co-0)')  smoothly  match 
the  second-order  bound-wave  solutions  of  dispersive  finite  depth  theory  (Herbers  and 
Burton,  1996). 

The  model  is  initialized  with  a  spectrum  E(co)  specified  at  the  offshore  boundary 
of  the  model  domain  (e.g.,  from  nearby  measurements  or  a  regional  model  prediction)  and 
the  corresponding  second-order  finite  depth  theory  expression  for  the  bispectrum  B(o)',o)- 

CO'). 
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flCo/.G)  -co')  =  2[D(g)/,g)-g)/)£((o/)£((«>-g>/) 


+  DCw^-q^Co^Co)  +  DCca-o^-w^w-w^Co))] 


(3) 


The  coupling  coefficient  D  is  given  by 


D(g)  ,w  )  = 

v   r  r 


(a)     +  G)  ) 

v    1  r 


#(k    +  K2)tanh[(Ki  +  k2)A]  -  (o>i  +  o>2) 


x   < 


G) 


0)  £K  K 

12  6     1    2 
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g  0)  0)         2(0)    +  o)  ) 

6  12  v     1  2y 


1 


K 


2  2 

0)  cosh  (k  h)      0)  cosh  (k  h) 

1  v    1    '  2  v    2    y 


/  J 


2  2 

0)        +G)G)+G)  £KK 

1  12  2  °    1    2 


(4) 


2g 


2o>  0) 

1      2 


where  /*  is  the  local  water  depth,  and  the  frequencies  0)j,  o)2  and  wavenumbers  Kj,  k2  (k, 
<  0  for  o)j  <  0)  of  the  interacting  primary  wave  components  obey  the  linear  theory 
dispersion  relation  0)t   =  g^ann^/i)  (Hasselmann,  1962;  Hasselmann  et  al.,1963). 
Using  the  symmetry  relations  (Hasselmann  et  al.,  1963): 

£(o)  =  £(-co) 
5(g)  ,co  -co  )  =  B(o)  -co  ,0) )  =5  *(-(o  ,G)  -co) 
=  B(o)  ,-o))  =  £(<•)-&)  ,-o>) 


where  *  indicates  the  complex  conjugate,  the  integral  term  in  equation  (1)  can  be 
expressed  as  the  sum  of  two  integrals  over  positive  frequencies 


0) 


jIM{B(co',co-co')}d^  =   flM{B((D/,(o-uy/)}d(o'  -  2r/M{fl(G)/,G>)}rfG>/  (5) 

-  «.  0  0 

that  represent  the  energy  transfers  to  frequency  co  resulting  from  the  sum  interaction  of 
frequencies   co',   co-co',  and  the  difference  interaction  of  frequencies      co+co',   co', 
respectively.  Hence  integrations  of  the  spectrum  and  bispectrum  evolution  equations  (1), 
(2)  can  be  restricted  to  positive  frequencies  (co',  co-co',  co  >  0). 
The  spectrum  and  bispectrum  are  discretized: 

0)    =  nAo)  for    n  =  1,  2,  ..,  N 

n 

E    =  E(co  )         for    n  =  1,  2,  ..,  AT 

n  n 

5      =  fl      +  *7     =  5(0)  ,0)  )    for    n  =  1,  2,  ..,  N-l  and    m  =  1,  2,  ..,  JV-i 

ran  ntn  nm  n       m 

where  A  co  is  the  bandwidth  and  coN  the  highest  frequency  included  in  the  computations. 
With  these  definitions,  equations  (1)  and  (2)  reduce  to  a  linear  set  of  N2  ordinary 
differential  equations  which  can  be  written  in  the  general  form 

dY  ^_ 
-r  =  F(Y) 
ax 

where  the  elements  of  y  are  the  discretized  spectrum  (En)  and  bispectrum  (Rnm,  Inm) 
values  and  F(Y)  incorporates  the  corresponding  right-hand  side  of  equations  (1)  and  (2). 
This  system  of  equations  is  solved  using  the  Bulirsch-Stoer  method,  a  variant  of 
Richardson  extrapolation  to  the  limit  that  uses  adaptive  stepsize  control  (Press  et  al., 
1992;  Bulirsch  and  Stoer,  1966). 
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III.  SIMULATIONS 

To  investigate  the  general  model  characteristics  and  the  dependence  of  wave 
shoaling  evolution  on  the  nonlinearity,  spectral  shape  and  bottom  profile,  extensive 
numerical  simulations  were  carried  out  with  simple  incident  wave  spectra  of  the  general 
form: 


/ 


£(/) 


n 
E— exp 


n 


1-h 


(6) 


E  CO  =  E sech 

JP 


J   JP 

n 


f. 


\ 


(7) 


where  E1  and  E2  are  single-sided  spectra  (E(f)  =  4izE((x))  with/=  o)/2ti),  E  is  the  surface 
elevation  variance  <r|  >  and  the  parameter  n  defines  the  width  of  the  spectrum.  E1  is  a 
broad  spectrum  with  a  power  law  high-frequency  tail  characteristic  of  actively  generated 
wind  waves.  E2  has  a  narrow  symmetric  exponential  shape,  characteristic  of  remotely 
generated  swell.  All  model  simulations  were  initialized  in  a  depth  A  =  6m  with  a  spectral 
peak  frequency^  =  0.07  Hz.  The  corresponding  wave  number  Kp  =  0.058  m"  yields  a 
representative  value  of  the  dispersion  parameter  Kh  =  0.35  at  the  offshore  boundary. 
Example  simulations  of  the  shoaling  of  a  broad  sea  spectrum  (E2  with  n  =  5,  the  Pierson- 
Moskowitz  [1964]  spectral  shape)  and  a  narrow  swell  spectrum  (E2  with  n  =  20,  the  full 
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width  at  half  maximum  power  is  0.009  Hz)  are  shown  in  Figures  1-5  for  different  bottom 
profiles  and  incident  significant  wave  heights  (Hs  =  4El/2).  The  initial  bispectrum  was 
obtained  by  substituting  the  initial  spectrum  in  the  finite  depth  theory  relation  (equation 
3). 

The  shoaling  evolution  of  a  narrow  swell  spectrum  with  significant  wave  heights 
of  0.05  and  0.5  m  on  plane  beaches  with  slopes  of  1:30  and  1:300  is  shown  in  Figure  1. 
All  four  simulations  show  the  familiar  growth  of  harmonic  peaks  at  frequencies  2fp,  3fp, 
••-  and  an  infragravity  peak  at  about  0.01  Hz.  Even  for  the  small  Hs  =  0.05  m  waves, 
harmonic  spectral  levels  are  significant  (up  to  10  %  of  the  primary  peak  level)  in  1.5  m 
depth.  Although  the  nonlinearity  remains  weak  ([2E\11  /  h,  a  representative  value  of  a/h, 
varies  between  0.003  and  0.017),  the  Ursell  number  (Ur  -  [2E]l/2 /  Kp2h3)  increases  from 
0.024  in  6  m  depth  to  a  relatively  large  value  of  0.58  in  1.5  m  depth.  As  expected,  the 
shoaling  evolution  is  much  stronger  for  the  larger  Hs  =  0.5  m  waves  with  harmonic 
spectral  levels  that  are  comparable  to  the  primary  peak  levels  in  1.5  m  depth.  In  these 
simulations  a/h  increases  from  0.03  in  6  m  depth  to  0.17  in  1.5  m  depth,  and  Ur  increases 
from  0.24  to  5.8.  In  both  the  Hs  =  0.05  and  0.5  m  cases,  stronger  growth  of  harmonic  and 
infragravity  peaks  is  predicted  on  a  gentle  1:300  slope  than  on  a  steep  1:30  slope. 
Eventually  (Figure  1  d,  f)  nonlinear  energy  transfers  fill  the  valleys  between  harmonic 
peaks  and  the  spectrum  flattens,  similar  to  simulations  with  a  deterministic  Boussinesq 
model  reported  by  Elgar  et  al.  (1990b). 

The  energy  transfers  and  their  dependence  on  the  bottom  slope  are  further 
illustrated  in  Figure  2  with  normalized  bispectrum  predictions  in  2  m  depth.  The 
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normalized  bispectrum: 
W,  /,)  -  =  (8) 


9 

with  B(fj,f2)  =  8tt  B(o)j,  o)2),  defined  analogous  to  £(/),  indicates  a  relative  measure  of 
phase  coupling  between  frequencies  fh  f2  and  f1  +  f2  (Herbers  et  al.,  1992).  Positive 
values  of  the  imaginary  part  of  the  bispectrum  indicate  energy  transfers  to  higher 
frequencies  through  sum  interactions,  whereas  negative  values  indicate  energy  transfers 
to  lower  frequencies  through  difference  interactions  (equations  1,5).  All  four  simulations 
show  strong  coupling  at  (/},  f2)  =  (0.07,  0.07)  Hz  associated  with  the  fp,  f ,  2fp  sum 
interaction  and  at  (0.07,  0.14)  Hz  (ihefp,  2fp,  3fp  sum  interaction).  The  larger  wave  and 
gentle  bottom  slope  simulations  also  show  coupling  to  higher  harmonics  (e.g.,  the  (0.14, 
0.14),  (0.21,  0.07),  and  (0.21,  0.14)  Hz  peaks).  Whereas  the  imaginary  part  of  b  is  small 
on  the  gentle  1:300  slope  (i.e.,  peaked  but  nearly  symmetric  wave  shapes,  Elgar  and 
Guza,  1985)  the  real  part  of  b  is  small  on  the  steep  1:30  slope  (i.e.,  pitched  forward  wave 
shapes).  The  energy  transfer  rate  (proportional  to  the  imaginary  part  of  the  bispectrum, 
equation  1)  is  smaller  on  the  gentle  slope  than  on  the  steep  slope  but  the  longer 
interaction  distances  cause  stronger  cumulative  spectral  evolution   (Figure  1). 

Simulation  results  of  the  shoaling  of  a  broad  spectrum  with  the  same  initial 
significant  wave  heights  (0.05  and  0.5  m)  and  beach  slopes  (1:30  and  1:300)  are  shown 
in  Figure  3.  The  spectral  evolution  is  much  weaker  than  in  the  narrow  spectra  simulations 
because  the  principal  effect  of  triad  interactions  is  to  spread  energy  to  frequencies  where 
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spectral  levels  are  relatively  low.  In  the  Hs  =  0.05  m  simulations  (Figure  3  a,  c,  e)  the 
nearly  uniform  increase  in  spectral  levels  at  frequencies  ^  0.05  Hz  is  a  linear  shoaling 

(decrease  in  wave  group  speed)  effect  ( the  first  term  on  the  right-hand  side  of  equation 
1).  The  1:300  slope  simulation  shows  slightly  larger  growth  of  spectral  levels  above  about 
2fp  that  is  the  result  of  sum  interactions.  The  larger  wave  (Hs  =  0.5  m)  simulations 
(Figure  3  b,d,f)  show  the  expected  stronger  nonlinear  evolution.  Although  the  spectrum 
remains  featureless,  nonlinear  interactions  cause  a  flattening  to  a  nearly  white  spectrum 
in  1.5  m  depth,  similar  to  the  narrow  spectrum  simulations  (Figure  If).  The  dependence 
on  bottom  slope  is  qualitatively  similar  to  the  narrow  spectrum  results  with  larger 
cumulative  energy  transfers  on  a  gently  sloping  beach  (compare  the  solid  and  dashed 
curves  in  Figure  3). 

The  shoaling  evolution  of  a  broad  spectrum  of  waves  (initial  Hs  =  0.5  m)  over 
three  different  bottom  profiles  is  compared  in  Figure  4.  All  three  profiles  start  with  a 
gently  sloping  (1:300)  section  from  6  to  3  m  depth  to  let  the  waves  evolve  to  a  shallow 
water  regime  with  significant  nonlinear  energy  transfers.  From  3  m  depth  shoreward,  the 
waves  either  continue  to  shoal  on  a  1:300  slope  to  1  m  depth  ("beach"  case),  propagate 
the  same  600  m  distance  in  (constant)  3  m  depth  ("flat"  case),  or  unshoal  over  a  -1:200 
section  back  to  6  m  depth  ("bar"  case).  In  contrast  to  the  beach  case,  the  spectral 
evolution  between  x  =  900  and  x  =  1500  m  is  weak  over  the  flat  bottom,  and  on  the  back 
side  of  the  bar  high  frequency  spectral  levels  are  reduced  to  approximately  the  initial 
levels  in  6  m  depth.  At  infragravity  frequencies  spectral  levels  continue  to  increase  on  all 
three  profiles,  but  the  growth  is  strongest  on  the  beach  profile  and  weakest  on  the  bar 
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profile  (Figure  4). 

Bispectra  after  150  m  of  evolution  on  the  three  different  profiles  are  compared  in 
Figure  5.  The  predicted  small  imaginary  part  of  the  bispectrum  remains  positive  on  the 
beach  profile,  allowing  for  continued  energy  transfers  to  high  frequencies  (equations  1, 
5).  On  the  flat  bottom  where  the  spectral  evolution  is  weak,  the  imaginary  part  of  the 
bispectrum  shows  small  alternating  positive  and  negative  peaks.  On  the  negative  slope 
section  of  the  bar  profile,  positive  values  of  the  imaginary  part  of  the  bispectrum  (i.e., 
energy  transfers  to  higher  frequencies)  evolve  to  negative  values  (i.e.,  energy  transfers  to 
lower  frequencies)  over  a  wide  range  of  frequencies,  causing  a  reversal  in  nonlinear 
energy  transfers  towards  lower  frequencies  (Figure  4). 

In  simulations  of  waves  propagating  over  a  flat  bottom  and  into  deeper  water 
(Figure  4),  small  undulations  are  noted  in  the  spectra  that  grow  with  distance  and 
eventually  develop  into  instabilities.  A  sensitivity  analysis  of  the  numerical  solutions  to 
variations  in  frequency  bandwidth,  the  error  tolerance  of  the  numerical  integration 
routine,  the  maximum  frequency,  and  using  different  extrapolation  techniques  (i.e., 
polynomial  and  rational  extrapolation),  yielded  identical  features  in  all  calculations. 
Energy  was  also  conserved  in  these  simulations  to  a  high  degree  of  accuracy.  These 
numerical  tests  indicate  that  the  predicted  growing  undulations  in  the  spectrum  are  true 
features  of  the  spectral  and  bispectral  evolution  equations  and  not  caused  by  numerical 
truncation  errors.  However,  the  Boussinesq  equations,  truncated  at  second-order  in 
nonlinearity  and  thus  valid  only  over  0(a/h)~  wavelengths  (Freilich  and  Guza,  1984),  do 
not  describe  accurately  the  longterm  evolution  of  these  moderately  energetic  waves. 
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Hence,  the  undulations  in  the  spectrum  may  not  be  physically  real  but  possibly  result 
from  the  breakdown  of  the  weakly  nonlinear  approximation  used  in  the  present  model. 
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IV.  COMPARISONS  TO  FIELD  OBSERVATIONS 

Extensive  field  observations  of  wave  shoaling  were  obtained  during  the  fall  of 
1994  in  the  DUCK94  nearshore  field  experiment  (Elgar  et  al.,  1996).  A  cross-shore 
transect  of  15  SPUV  systems,  each  consisting  of  a  co-located  pressure  transducer, 
bidirectional  electromagnetic  current  meter,  and  sonar  altimeter  were  deployed  on  a  sandy, 
barred  beach  near  Duck,  North  Carolina.  The  350  m  long  transect  extended  from  the 
shoreline  to  about  6  meters  depth  (Figure  6).  The  sample  frequency  of  all  instruments  was 
2  Hz.  Sea  surface  elevation  spectra  with  approximately  120  degrees  of  freedom  were 
estimated  from  three-hour-long  pressure  records  using  a  linear  theory  depth  correction. 

The  present  analysis  of  four  case  studies  is  focused  on  benign  wave  conditions 
(incident  wave  significant  heights  ranged  from  0.4  -  0.8  m)  when  the  surf  zone  was 
confined  to  the  beach  face  at  the  shoreward  end  of  the  transect.  These  observations  span 
a  2  week  period  in  September  with  small  bathymetric  changes.  Differences  between  the 
depth  profiles  of  the  case  studies  (Figure  6)  are  primarily  due  to  tidal  sea  level 
fluctuations.  The  beach  profile  is  characterized  by  a  sandbar  located  about  120-140  m 
from  the  shoreline  and  submerged  approximately  2.2-2.5  m  below  the  mean  sea  surface. 
The  bottom  slope  is  approximately  1:80  seaward  of  the  sandbar.  Shoreward  of  the 
sandbar,  the  profile  deepens  slightly  (20-40  cm)  into  a  relatively  flat,  80  m  wide  trough 
that  extends  to  the  steep  (1:10)  beach  face.  The  beach  profile  used  in  the  Boussinesq 
model  computations  was  obtained  through  linear  interpolation  of  the  depth  estimates 
extracted  from  the  sonar  altimeters  (Gallagher  et  al.,  1997). 
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Stochastic  Boussinesq  model  predictions  of  wave  spectrum  evolution  in  the  four 
case  studies  are  compared  to  the  observed  spectrum  evolution  and  deterministic 
Boussinesq  model  predictions  (provided  by  Dr.  Steve  Elgar)  in  Figures  7,  8,  12  and  15. 
The  stochastic  model  predictions  for  15,  21  and  24  September  were  initialized  with  the 
spectrum  measured  at  the  furthest  offshore  pressure  sensor,  pi 9.  The  10  September  case 
was  initialized  with  pressure  sensor  pi 8  because  sensor  pl9  malfunctioned  during  this 
run.  The  deterministic  model  predictions  were  initialized  with  the  measured  pressure  time 
series  (see  Elgar  et  al.  1996  for  further  details)  at  the  same  locations. 

In  all  four  cases  (and  other  case  studies  not  shown)  predictions  of  both  models 
are  in  excellent  agreement  with  the  observed  wave  shoaling  evolution.  The  narrow  swell 
spectrum  cases  (10  and  15  September)  show  the  familiar  amplification  of  harmonic  peaks. 
On  10  September  the  incident  wave  spectrum  is  dominated  by  a  narrow  swell  peak  (peak 
frequency^  ~  0.075  Hz),  with  a  broader,  but  relatively  small  sea  peak  at  0.12  Hz  (Figure 
7).  Energy  is  transferred  from  the  swell  peak  frequency  fp  to  higher  frequencies  through 
sum  triad  interactions  resulting  in  distinct  harmonic  peaks  at  2fp  (0.15  Hz;  driven  by^, 
fp  interactions),  3fp  (0.23  Hz;  fp,  2fp  interactions)  and  4fp  (0.3  Hz;  fp,  3fp  and  2fp,  2fp 
interactions).  Close  to  shore  the  small  0.12  Hz  incident  sea  peak  is  completely  submerged 
in  the  2fp  swell  harmonic  (Figure  7f). 

On  September  15  the  incident  wave  spectrum  is  distinctly  bimodal  with  a  narrow 
swell  peak  (fp  ~  0.06  Hz)  and  a  slightly  smaller  sea  peak  at  twice  the  swell  frequency  (2fp 
=  0.12  Hz)  (Figure  8).  As  these  wave  systems  propagate  over  the  shallow  sandbar,  large 
nonlinear  energy  transfers  in  sum  interactions  yield  clearly  distinguishable  harmonic  peaks 
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at  3fp,  Afp  and  5fp  (compare  Figure  8a  and  8c).  The  cross-shore  evolution  of  spectral  levels 
at  frequencies  fp,  2fp,  3fp,  4fp  and  5fp  is  shown  in  Figure  9.  Partial  reflection  of  the 
dominant  0.06  Hz  swells  from  the  beach  is  evident  in  the  large  cross-shore  energy 
variations  (i.e.,  standing  wave  patterns)  observed  close  to  shore.  Good  agreement  between 
the  observed  and  predicted  growth  of  higher-frequency  harmonics  indicates  that  nonlinear 
energy  transfers  are  insensitive  to  weak  reflections  from  shore  (Elgar  et  al.,  1996). 

Both  the  observed  and  predicted  bispectra  on  September  15  show  the  expected 
shoaling  transition  from  real  values  (i.e.,  skewed  but  symmetric  wave  profiles)  (Figure 
10a  and  10b)  to  imaginary  values  (i.e.,  pitched-forward  wave  profiles)  (Figure  10c  and 
lOd).  Although  the  observed  and  predicted  bispectral  levels  generally  agree  within  the 
considerable  uncertainty  in  estimates  extracted  from  short  field  data  records  (records 
longer  than  3  hours  could  not  be  used  owing  to  tide-induced  depth  changes),  they  differ 
in  detail  at  the  shallower  sites  (Figure  1 1).  In  frequency  pairs  involving  the  0.06  Hz  swell 
peak,  the  observed  bispectrum  shows  a  dramatic  shift  from  real  to  imaginary  values 
between  sensors  p3  and  p4  (separated  by  only  15  m)  that  is  absent  in  the  model  results. 
This  biphase  shift  is  caused  by  the  partial  reflection  of  the  0.06  Hz  swell  from  shore 
(Figure  9)  that  is  not  incorporated  in  the  model  predictions  .  Midway  between  nodes  and 
antinodes  (e.g.,  sensor  p4)  the  reflected  components  are  90°  out  of  phase  with  the  incident 
components,  causing  large  biphase  shifts  in  triads  involving  the  standing  wave  component. 

The  September  24  spectrum,  is  broader  and  more  energetic,  with  a  peak  frequency, 
fp  ~  0.1  Hz  (Figure  12).  Sum  interactions  transfer  energy  to  a  broad  range  of  higher 
frequencies,  causing  a  broadening  of  the  spectrum  rather  than  the  development  of  distinct 
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harmonic  peaks  observed  on  September  10  and  15.  The  observed  and  predicted  spectral 
levels  at  3fp  decrease  sharply  between  the  bar  crest  (x  =  240  m)  and  the  slightly  deeper 
trough  (x  =  300  m),  whereas  energy  levels  at  2fp  continue  to  increase  (Figure  13).  These 
results  suggest  that  energy  is  transferred  back  to  lower  frequencies  as  waves  travel  over 
the  sandbar  into  deeper  water  similar  to  the  simulation  results  (Figures  4,  5).  Observed 
and  predicted  bispectra  (Figure  14)  show  a  clear  transition  from  positive  imaginary  parts 
seaward  of  the  bar  crest  to  negative  imaginary  parts  shoreward  of  the  bar  crest  at  high 
frequencies  that  is  consistent  with  this  reversal  in  the  nonlinear  energy  transfer  (equations 
1,5). 

In  contrast  to  the  September  10,  15  and  24  case  studies,  the  shoaling  evolution  of 
the  broad,  featureless  spectrum  observed  on  September  21  (Figure  15)  is  weak  even 
though  the  nonlinearity  is  comparatively  strong  in  this  case  (the  incident  wave  significant 
height  is  0.8  m).  Sum  and  difference  interactions  cause  significant  energy  transfers,  but 
do  not  strongly  affect  the  spectrum  because  the  interactions,  spread  over  a  wide  frequency 
range,  tend  to  cancel  out  in  a  broad  spectrum.  Predicted  bispectral  levels  (not  shown)  are 
low,  consistent  with  the  observations. 

Discrepancies  between  the  stochastic  model  predictions  and  the  observations,  and 
between  the  deterministic  model  predictions  and  the  observations  are  roughly  comparable 
(observed  and  predicted  spectral  levels  generally  agree  within  about  a  factor  of  four)  but 
differ  in  detail.  Initially  the  deterministic  model  tends  to  overpredict  energy  transfers  to 
higher  frequencies  whereas  the  stochastic  model  predictions  are  in  close  agreement  with 
the  observed  spectra  (e.g.,  Figures  7a-c,  8a-b).  During  the  later  stages  of  shoaling 
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evolution,  the  stochastic  model  tends  to  overpredict  high-frequency  spectral  levels, 
whereas  the  deterministic  model  predictions  are  close  to,  or  in  some  cases,  underpredict 
(e.g.,  Figure  12  c,  d)  the  observed  spectral  levels.  Some  of  these  differences  may  be  the 
result  of  the  different  formulation  of  dispersion  effects  in  the  two  models.  The  stochastic 
model  uses  Peregrine's  (1967)  "consistent"  approximation,  whereas  the  deterministic 
model  uses  Freilich  and  Guza's  (1984)  "dispersive"  approximation.  Comparisons  of 
deterministic  Boussinesq  models  to  field  data  reported  by  Freilich  and  Guza  (1984)  show 
similar  trends  of  overprediction  (for  the  "consistent"  approximation)  and  underprediction 
("dispersive"  approximation)  of  high  frequency  spectral  levels.  Other  possible  explanations 
for  small  differences  between  the  deterministic  and  stochastic  model  predictions  are  the 
different  way  the  models  initialize  third-order  statistics  (the  stochastic  model  uses  second- 
order  finite  depth  theory  whereas  the  deterministic  model  uses  measured  time  series)  and 
the  statistical  closure  of  the  stochastic  model.  Additionally,  dissipation  (neglected  in  both 
models)  and  higher-order  nonlinear  effects  likely  contribute  significant  errors  in  the 
predictions  close  to  shore. 

The  predicted  shoaling  amplification  of  low  frequency  (<0.06  Hz)  spectral  levels 
is  generally  in  reasonable  agreement  with  the  observations,  even  though  the  model  is 
obviously  inadequate  at  infragravity  frequencies.  It  is  well  known  that  reflection  from 
shore  (e.g.,  Elgar  et  al.,  1994)  and  refractive  trapping  in  deeper  water  (e.g.,  Huntley  et  al., 
1981;  Herbers  et  al.,  1995  a,  b)  are  important  at  infragravity  frequencies,  and  these  effects 
are  not  incorporated  in  the  models.  Furthermore,  the  energy  transfers  to  infragravity 
frequencies  are  sensitive  to  directional  spreading  effects  that  are  neglected  here  (Herbers 
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et  al.,  1995;  Herbers  and  Burton,  1996).  Nevertheless,  the  roughly  comparable  observed 
and  predicted  infragravity  energy  levels  suggest  that  nonlinear  triad  interactions  is  a 
plausible  mechanism  for  the  transfers  of  energy  to  infragravity  frequencies  in  shallow 
water. 
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V.  SUMMARY  AND  CONCLUSIONS 

A  stochastic  model  for  the  shoaling  of  waves  on  a  beach  with  straight  and  parallel 
depth  contours  is  presented,  based  on  a  third-order  closure  of  Boussinesq  equations 
(Herbers  and  Burton,  1996).  The  model  includes  nonlinear  triad  interactions  in  which  two 
primary  wave  components  with  frequencies  (Oj  and  co2,  excite  a  secondary  wave 
component  with  the  sum  ((Oj+O)^  or  difference  (o)j-0)2)  frequency.  Neglecting  directional 
spreading  effects,  a  coupled  set  of  evolution  equations  for  the  wave  spectrum  and 
bispectrum  is  solved  with  standard  numerical  integration  techniques.  The  model  is 
numerically  efficient  and  requires  only  an  estimate  of  the  incident  wave  spectrum  for 
initialization  that  is  often  readily  available  from  a  nearby  buoy  or  a  regional  wave  model 
prediction.  The  bispectrum  is  initialized  with  a  local  prediction  based  on  second-order 
finite  depth  theory. 

Extensive  numerical  simulations  were  performed  to  examine  the  model 
characteristics  and  the  dependence  of  wave  shoaling  evolution  on  nonlinearity,  spectral 
shape  and  bottom  profile.  In  simulations  with  strong  nonlinearity,  both  narrow  and  broad 
spectra  tend  to  evolve  to  a  flat  featureless  spectrum  (Figures  If,  3f).  Simulations  of 
narrow  spectra  with  peak  frequency^  show  the  familiar  growth  of  harmonic  peaks  at  2fp, 
3fp,  ••-.  In  simulations  with  broad  spectra  comparable  energy  transfers  to  higher 
frequencies  occur,  but  since  the  interactions  are  spread  over  a  wide  frequency  range,  the 
spectra  remain  featureless  (Figure  3).  On  gently  sloping  beaches  predicted  nonlinear 
energy  transfer  rates  are  weaker  than  on  steep  beaches  but  cause  stronger  cumulative 

23 


spectral  evolution  (Figures  1,3).  On  steep  slopes,  predicted  bispectra  have  relatively  large 
imaginary  parts  characteristic  of  pitched  forward  wave  shapes,  whereas  the  predominantly 
real  bispectral  values  predicted  on  gentle  slopes  indicate  symmetric  wave  profiles  .  These 
characteristics  are  qualitatively  consistent  with  wave  shape  evolution  observed  prior  to 
breaking  on  natural  beaches.  Simulations  of  waves  propagating  over  a  bar  into  deeper 
water  show  decreasing  high-frequency  spectral  levels  (Figure  4).  The  predicted  bispectra 
indicate  a  reversal  in  nonlinear  energy  transfers  with  difference  triad  interactions 
transferring  energy  back  towards  lower  frequencies. 

Stochastic  and  deterministic  (Freilich  and  Guza,  1984)  Boussinesq  model 
predictions  were  compared  to  extensive  field  observations  of  wave  shoaling  on  a  natural 
barred  beach.  Although  predictions  of  the  two  models  differ  in  detail,  the  overall 
agreement  with  the  observed  wave  spectrum  evolution  is  comparable.  Both  models  predict 
accurately  the  nonlinear  transfer  of  energy  to  higher  frequencies  for  a  range  of  incident 
wave  conditions  (Figures  7,  8,  12,  15).  These  results  are  similar  to  earlier  studies  using 
deterministic  Boussinesq  models  on  nearly  plane  California  beaches  (Freilich  and  Guza, 
1984;  Elgar  and  Guza,  1985).  Although  spectral  levels  at  high  frequencies  generally 
increase  as  waves  propagate  shoreward  owing  to  sum  triad  interactions,  in  one  case  a 
decrease  in  high-frequency  spectral  levels  was  observed  shoreward  of  the  sandbar, 
consistent  with  difference  interactions  predicted  by  both  models  (Figures  12,  13).  These 
observations  support  the  simulation  result  that  energy  can  be  transferred  back  to  incident 
wave  frequencies  in  regions  of  gradually  increasing  depth. 

In  conclusion,  a  numerically  efficient  stochastic  Boussinesq  model  was  developed 
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that  predicts  the  evolution  of  wave  spectra  and  bispectra  on  beaches.  Comparisons  of 
model  predictions  to  extensive  field  observations  from  a  natural  beach  show  excellent 
agreement.  In  the  future,  the  stochastic  wave  shoaling  model  will  be  extended  to 
directionally  spread  waves  based  on  similar  evolution  equations  for  the  frequency- 
alongshore  wavenumber  spectrum  and  bispectrum  (Herbers  and  Burton,  1996). 
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APPENDIX 

Figure  1.  Numerical  simulations  of  the  shoaling  evolution  of  a  narrow  spectrum 
of  waves  (equation  7  with  n  =  20)  over  a  plane  beach.  The  model  is  initialized  in  6  m 
depth.  Results  are  shown  in  4  m  (upper  panels),  2  m  (middle  panels),  and  1.5  m  (lower 
panels)  depth  for  incident  wave  significant  heights  of  0.05  m  (left  panels)  and  0.5  m 
(right  panels)  and  beach  slopes  of  1:300  (solid  lines)  and  1:30  (dashed  lines).  The  initial 
spectrum  is  indicated  in  each  panel  with  a  dotted  line. 

Figure  2.  Normalized  bispectra  b(fj,f2)  (equation  8,  units  Hz" '  )  predicted  in  2 
m  depth  in  the  simulations  described  in  Figure  1.  The  real  and  imaginary  parts  of  b  are 
shown  in  the  lower  and  left  quadrants,  respectively.  Contour  levels  are:  ±  1,  3,  5. 

Figure  3.  Numerical  simulation  of  the  shoaling  evolution  of  a  broad  spectrum  of 
waves  (equation  6  with  n  =  5)  over  a  plane  beach  (same  format  as  Figure  1). 

Figure  4.  Simulation  of  the  evolution  of  a  broad  spectrum  (equation  6  with  n  = 
5),  with  an  initial  Hs  =  0.5  m  in  6  m  depth,  over  3  different  bottom  profiles:  a  beach  (left 
panels),  a  flat  section  (center  panels),  and  a  bar  (right  panels).  Predictions  are  shown  at 
x  =  900,  1050,  and  1200  m  (indicated  by  asterisks  in  the  upper  panels).  The  dotted  line 
indicates  the  initial  spectrum  at  x  =  0. 

Figure  5.  Normalized  bispectra  predicted  in  the  simulations  described  in  Figure 
4.  (a)  x  =  900  m  (same  for  all  3  profiles),  (b)  x  =  1050  m  on  the  beach  profile,  (c)  x  = 
1050  m  on  the  flat  section,  (d)  x  =  1050  m  on  the  bar  profile.  The  format  of  the  panels 
is  the  same  as  in  Figure  2. 
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Figure  6.  Depth  profiles  and  sensor  locations  of  the  4  field  data  case  studies. 

Figure  7.  Comparison  of  observed  (solid  line)  and  predicted  (asterisks  =  stochastic 
model,  circles  =  deterministic  model)  spectra  on  September  10  at  6  instrument  locations. 
The  initial  spectrum  (Hs  =  0.5  m)  is  indicated  in  each  panel  with  a  dashed  line. 

Figure  8.  Comparison  of  observed  and  predicted  spectra  on  September  15  (Hs  = 
0.4  m)  (same  format  as  Figure  7). 

Figure  9.  Cross-shore  evolution  of  the  spectral  levels  at  the  peak  frequency  and 
the  first  four  harmonics  on  September  15.  The  solid  lines  are  the  observed  levels,  the 
dotted  lines  are  the  detenninistic  model  predictions,  and  the  dashed  lines  are  the  stochastic 
model  predictions. 

Figure  10.  Bispectra  observed  (right  panels)  and  predicted  (left  panels,  from  the 
stochastic  model)  offshore  of  the  bar  (sensor  pi 7,  upper  panels)  and  on  the  bar  (sensor 
pl4,  lower  panels)  on  September  15  .  The  format  of  the  panels  is  the  same  as  in  Figure 
2. 

Figure  11.  Bispectra  observed  (right  panels)  and  predicted  (left  panels,  from  the 
stochastic  model)  on  September  15  at  sensors  p4  (upper  panels)  and  p3  (lower  panels), 
both  located  on  the  beachface  close  to  shore.  The  format  of  the  panels  is  the  same  as  in 
Figure  2. 

Figure  12.  Comparison  of  observed  and  predicted  spectra  on  September  24  (Hs 
=  0.8  m)  (same  format  as  Figure  7). 

Figure  13.  Cross-shore  evolution  of  the  spectral  levels  at  the  peak  frequency  and 
the  first  two  harmonics  on  September  24  (same  format  as  Figure  9). 
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Figure  14.  Bispectra  observed  (right  panels)  and  predicted  (left  panels,  from  the 
stochastic  model)  on  the  bar  (sensor  14,  upper  panels)  and  inshore  of  the  bar  (sensor  p23, 
lower  panels)  on  September  24  .  The  format  of  the  panels  is  the  same  as  in  Figure  2. 

Figure  15.  Comparison  of  observed  and  predicted  spectra  on  September  2 1  (Hs 
=  0.8  m)  (same  format  as  Figure  7). 
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